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ABSTRACT 

We use time-evolutions of tlic linear perturbation equations to study the oscillations 
of rapidly rotating neutrons stars. Our models account for the buoyancy due to com- 
position gradients and we study, for the first time, the nature of the resultant g-modes 
in a fast spinning star. We provide detailed comparisons of non-stratified and strati- 
fied models. This leads to an improved understanding of the relationship between the 
inertial modes of a non-stratified star and the g-modes of a stratified system. In partic- 
ular, we demonstrate that each g-mode becomes rotation-dominated, i.e. approaches 
a particular inertial mode, as the rotation rate of the star is increased. We also dis- 
cuss issues relating to the gravitational-wave driven instability of the various classes 
of oscillation modes. 
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1 INTRODUCTION 



The oscillations of rotating neutron stars are of great interest in astrophysics. Recent efforts to understand the various pulsation 
modes of compact stars have to a large extent been motivated by gravitational-wave astronomy ( Andersson fc Kokkotas , IQoil . 
1998l V Different modes of oscillation depend on different pieces of internal physics and one may hope to use o bservations to 



probe, for example, the composition of the high-densi tv region |Ferrari fc Gualtieri 20081 : Benhar et al. 2007 ). Work in this 



area, which has a long history (for a review see e.g., Kokkotas fc Schmidtl 1993), gained further momentum recently with 



the relatively successful matching between ob served quasiperiodic osci l lations in the tails of magnetar flares and calculated 
torsio nal oscillations of the neutron star crust (Duncan 19981 : Piro 2005 : Watts fc Strohmavei 2007 : Samuelsson fc AnderssonI 
20071). The excitement following these, the first ever, observations of likely neutron star vibrations is significant. It is clear that 
our models need to be improved significantly if we expect to carry out a quantitative astero-seismology programme for neutron 
stars, but the prospects for improvements look good. After all, we already have a r easonable understanding of the dynamics 
of the crust region ( see Carter fc Samuelsson 1 20061 : Samuelsson fc Andersson 2007) as well as the superfiuid components in 



the core of the star (Andersson fc Comei 



200' 



The present work is motivated by the need to consider more realistic models of rapidly rotating stars. Our emphasis 
will be on the role of the internal stratification a ssociated with composition variations. This leads to the presence of the 
so-called g-modes ( Reisenegger fc Goldreicbl 1992), and we want to investigate how these modes are affected by fast rota- 
tion (the analogous problem fo r therm al ocean g-modes of rapidly rotating neutron stars has already been considered, see 
Bildsten. Ushomirskv fc Cutler ( 19961 )). In order to determine the rotational effects on the oscillation spectrum we study 
linear perturbations of axisymmetric stellar configurations, where the Coriolis force and the centrifugal flattening of the star 
are included in the background. In this way, we go beyond the so-called slow rotation approximation, where the rotation 
itself is treated perturbatively and the mode-frequencies of non-rotating models are determined by a perturbation expansion 
with respect to the angular velocity of the star Q.. We are particularly intereste d in the relationship between the g-modes 
of a stratified stellar model and the inertial modes ( Lockitch fc Friedman 19991 ). for which the Coriolis force is the main 
restoring agent. The inertial modes are important since they ma y be driven unstable by the emission of gravita tional waves. 
The strongest instability is associated with the so-called r-modes ( Andersson fc KokkotasI 200ll : Andersson 20031 ). A key issue 
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concerns the amplitude at which an unstable mode saturates. De tailed work has shown that an unstable r-mode saturates due 
to nonlinear couplings to other inertial modes I Arras et al. 2003l l . However, the mode saturation has only been considered for 
slowly rotating, non-stratified models and it is important to ask to what extent these results will change if the models are 
made more re alistic. In stratif i ed neutron stars, the g-modes ten d to become rotationally dominated beyond a cert ain rotation 
rate (see e.g., Yoshida fc Lee 2000l : Dintrans &: RieutordI 200oh and may then be called inertia-gravity modes ( Unno et al. 
19891 ). Effects due to the stratification are thus weakened, which indicates that the r-mode saturation estimates may still be 



reliable, at least for the fastest spinning systems. 

I mproving and extending a numerical code that has been used to study non-stratified stars IjJones. Andersson fc Stergioulas 
20021 ). we approach the problem via time-evolutions of the equations that govern the linear perturbations of a, potentially 



rapidly, rotating Newtonian star. This strategy has the advantage that one does not have to deal with the many rotationally 
coupled multipoles, that tend to make a slow-rotation calculation less tractable. On the other hand, a numerical evolution 
is not expected to have the precision of a frequency domain calculation. Neither should one expect it to yield the complete 
spectrum of modes. After all, the simulation results depend on the chosen initial data. In absence of a clear idea of the nature 
of the various modes one does n ot know how to excite specific oscillatio ns. One way around this problem is t o follo w the 
"recycling" strategy developed by Stergioulas. Apostolatos fc Font 1 2004h and Dimmelmeier. Stergioulas fc Font ( 20061 ). An- 
other, more pragmatic, approach is to simply study the modes that are excited by "generic" initial data. This is the attitude 
that we adopt here. 



2 ROTATING STELLAR MODELS 

We want to investigate the effect that composition variations have on the various oscillation modes of a rapidly rotating neutron 
star. To do this, we consider an equilibrium configuration representing a rapidly and uniformly rotating star containing three 
particle species. We assume that matter is composed of degenerate neutrons n, protons p and electrons e. For simplicity, we 
further assume that neutrons and protons are not superfluid. In essence, our model represents the conditions in the outer 
core of a very young neutron star that has not yet cooled below the temperature at which the baryons form superfluid 
condensates (in practice with a core temperature significantly above 10^ K). We also do not consider the various exotic phases 
of matter that may be present in the deep core of a realistic neutron star model, e.g. hyperons or deconfined quarks (see 
Lattimer fc Praka sh"2004'. '2007^. Accounting for the presence of such phases would not be too difficult, but at this stage our 



focus is on understanding the basic role of composition variations. We are already extending this work to account for the 
presence of superfluid components and expect to report on results in this direction soon. 

The matter composition is represented by the three number densities rix, where x is either n, p or e. Charge neutrality is 



provi ded by the Coulomb interaction, that efficiently lock together protons and electrons (e.g.. IVillain. Bonazzola fc Haensel 



20051 ). Thus we have n-p — ric and Vp = Vc, where Vx denotes the velocity of each fluid component. In this model, the pressure 



P is a function of two variables, e.g. the neutron and proton number densities or alternatively the total baryon number density 
nb = rbn + rip and the proton fraction Xp = rip/rib. In the unperturbed stellar model, it is assumed that the weak interaction 
processes: 

n — > p + e + i/e, p + e — > n + , (1) 

lead to /3-equilibrium. The proton fraction is then only a function of the baryon number density (or equivalently the baryon 
mass density p — m (rip + Un), where for simplicity the proton and neutron masses have been assumed equal). Then the stellar 
fluid can be described by a barotropic equation of state P = P (p). In this paper, we consider a simple polytropic model 

P = Kp^P , (2) 

where K is constant and 

dlogP 



- 13 



(3) 



dlogp 
is the adiabatic index. 

Once the equation of s tate is provided, rotating stellar conflgurations can be de termined by the H achisu self-consistent 
fleld method ( Hachisu 19861 ). For this purpose we use the numerical code described in Jones et al. ( 2002 ). This code generates 



axisymmetric equilibrium models by specifying the polytropic index and the ratio of the polar to equatorial axes Rp I Req • 
In Table [T] we give some details relating to the rotation rate, energies and masses of equilibrium configurations used in this 
paper. 



3 THE PERTURBATION PROBLEM 

In a perturbed configuration, an oscillating fluid element may not have time to reach /3-equilibrium with the neighbouring 
matter. In fact, for typical stellar oscillation frequencies the time scale of the weak interaction processes is much too long for 
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Table 1. This table displays the main quantities of our rotating equilibrium configurations. The stellar models are described by a 
polytropic equation of state with adiabatic index Fp = 2. In the first and second columns we show, respectively, the ratio of polar 
to equatorial axes and the angular velocity of the star. In the third column, the rotation rate is compared to the Kepler velocity Qj^ 
that represents the mass shedding limit. The ratio between the rotational kinetic energy and gravitational potential energy T/|iy| and 
the stellar mass are given in the fourth and fifth columns, respectively. All quantities are given in dimensionless units, where G is the 
gravitational constant, pc represents the central mass density and R^q is the equatorial radius. 



Rp/B.eg n/V^ O/Qk T/\W\ X 102 M/{p,Rlq) 



1.000 


0.000 


0.000 


0.000 


1.273 


0.996 


0.084 


0.116 


0.096 


1.270 


0.992 


0.119 


0.164 


0.192 


1.261 


0.983 


0.167 


0.230 


0.385 


1.248 


0.950 


0.287 


0.396 


1.170 


1.197 


0.933 


0.331 


0.456 


1.570 


1.171 


0.900 


0.403 


0.556 


2.380 


1.118 


0.850 


0.488 


0.673 


3.640 


1.038 


0.800 


0.556 


0.767 


4.930 


0.956 


0.750 


0.612 


0.844 


6.250 


0.869 


0.700 


0.658 


0.907 


7.570 


0.779 


0.650 


0.693 


0.956 


8.820 


0.684 


0.600 


0.717 


0.989 


9.860 


0.579 



the matter to equilibrate within one oscillation period (e.g.. iReisenegger fc Goldreichlll992f ). In the limit of slow reactions one 
can thus accurately assume that the composition of a perturbed fluid element is frozen. This condition can be imposed by 
letting the Lagrangian variation of the proton fraction vanish, i.e., 

Aa;p = (5a;p + C • Vxp = , (4) 

where ^ denotes the Lagrangian displacement of a fluid element. 

The perturbed fluid is no longer described by the one-parameter equation of state that determined the background model. 
The perturbed pressure now depends on both the proton fraction and the baryon density, 

5P = 5P{p,x^). (5) 

Given equation Q, the Lagrangian perturbation of pressure is given by 



where 



AP^Ff-Ap, 
P 



aiogP 



(6) 



(7) 



91ogp 

is the adiabatic index for a frozen composition. Alternatively, equation (|6]) can be written in terms of Eulerian perturbations: 

5p _ 1 5P 

The vector field A is defined as 



Vlnp VlnP= VlnP 



(8) 



(9) 

Its magnitude is usually referred to as the Schwarzschild discriminant, with an overall sign depending upon whether A points 
outward (positive discriminant) or inward (negative discriminant). Note that the background relation ([3]) has been used in 
the last equality of equation (|9]). 

For barotropic stellar oscillations a single equation of state describes both the background and the perturbed matter. The 
adiabatic indices F^ and Ff are then equal and the Schwarzshild discriminant vanishes, A = 0. In this case, the stellar structure 
does not sustain buoyancy restored modes so the composition g- modes are absent from the spectrum IjReisenegger fc Goldreich 
19921 ). In the nonbarotropic case, A 7^ 0, the sign of the Schwarzschild discriminant determines the convective stability of 



the motion. Stab le (unstable) g modes are present in stars with negative (positive) Schwarzschild discriminant (|TassouL,19"78l : 
Unno et al.lll989l '). 



In our numerical evolutions, we do not use the Lagrangian displacement ^ as a dynamical variable. Instead we evolve the 
perturbed proton fraction 5xp. The relevant dynamical equation follows if we rewrite equation ([8| using equation (|4]); 



5p 
P 



_1_5P 

r; p" 



p 



(10) 
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where we have defined 



We have also assumed that Xp = Xp (p). 



Sxp/i'-^ 



(11) 



The composition gradient can be directly related to the deviation of Ff from the background Tp. From equation ([5| the 
pressure perturbation can be expressed as 



5P_ 
T 



5p 



Sxp 



(12) 



where 



dlogP 



d log Xr- 



(13) 



Comparing equations (|10p and (|12p we obtain the following relation between the proton fraction and the adiabatic indices 



^ dlogp 



(14) 



It is now straightforward to extend the models considered bv I Jones et alj (|2002f ) to account for composition variations. 
Non-axisymmetric perturbations of rapidly, and uniformly, rotating stratified neutron stars can be described by a system of 
five equations. This system is formed by the perturbed Euler, mass conservation and frozen composition equations. In the 
frame of the rotating background and in terms of the dynamical variables 5P, 5xp ^'^'^ the fiux f = pv, the equations can be 
written 

VP, 



dti = -V(5P - 2 (fi X f) 



VlnP 



<5P+ 1 



Tf 



dtSP = 



P 1 
^Tf- V - f + - 
P P 

f ■ VP 



ft-' 



f ■ VP, 



(15) 
(16) 
(17) 



In writing down equation (|15|l we have made the Cowling approximation, i.e. neglected the forces due to the perturbed 
gravitational potential. Equation (|17p describes the conservation of the proton fraction in a fluid element moving with the 
fluid. It is equivalent to the time derivative of equation ll4ll . For V g = Ff, the system of equations (|15|l - (|17|l reduces to the 
barotropic perturbation equations considered in I Jones et al. I l|2002l ). In that case, Sxp is a superfluous dynamical variable. 

As is clear from the above equations, time evolution of the non-axisymmetric perturbation equations is a three-dimensional 
problem in space. However, by exploiting the axial symmetry of the stellar background and the fact that we only consider 
linear perturbations, all perturbe d variables can be expanded in terms of a set of basis functions (cos mi;/) , sin mi;/)), where m 
is the azimuthal harmonic index i Papaloizou fc Pringle 1980l ). For instance, in the case of the pressure we have 



5P{t,r,e,(t)) = ^ [5 P:^^ {t, r, e) COS m(t) + 5 P~{t,r,e) sin ttkI)] 



(18) 



The perturbation equations then decouple with respect to m and the problem becomes two-dimensional. For any m, the 
non-axisymmetric oscillations of a stratifled neutron star requires the solution of a system of ten partial differential equations 
for the ten variables (5P^ , f * , 5Xp) ■ 



4 BOUNDARY CONDITIONS 

The evolution equations need to be complemented by boundary conditions. At the stellar surface we require that the La- 
grangian perturbation of the pressure vanishes, i.e., 

AP = 5P + C- VP = 0. (19) 

From the Euler equation for the stationary background 

VP = -pfi X (n X r) - pV$ , (20) 

the condition (|19[) reduces to 

AP = 5P- p^ • (n X n X r + V*) = 0. (21) 

In a polytropic star, the mass density vanishes on the stellar surface. Therefore, if the Lagrangian displacement ^ assumes 
flnite values at the surface, according to equation (|21|l the pressure perturbation i5P behaves as p near the surface. In order 
to avoid divergent solutions in the numerical simulations, we impose a zero surface condition for the mass flux perturbation 

f = pv = . (22) 
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We consider initial conditions where the Lagrangian displacement is finite on the stellar surface. Therefore, equation (|22|) 
implies that along the numerical evolution also the quantity remains zero on the surface. As a result, from equation pip 
the Eulerian pressure perturbation satisfies a zero surface condition 

SP^O. (23) 

This condition ensures that the perturbed quantities remain regular at the surface. It is also straightforward to implement in 
the code. 

We focus on the non-axisymmetric oscillations for m ^ 2. At the origin, r = 0, the regular behaviour of equations (I15|l - (|17| l 
is then guaranteed by the conditions: 

5P = <5xp = f = . (24) 



5 OSCILLATION MODES OF ROTATING STARS 

Before we discuss the results of our evolutions, it is useful to digress on the expected nature of the pulsation modes of a 
spinning star. 

Perturbations of spherical (non-spinning) stars can be decomposed in two separate classes identified by their behaviour 
under the parity transformation. Polar perturbations, which are expressed in terms of the scalar spherical harmonics Yj™ and 
their gradients \7Yi,„, transform as ( — 1)' under parity inversion. Meanwhile, axial perturbations, which behave as r x VYim, 
change sign as (—1)'"''^ under parity. In a spherical star, the different multipoles do not couple. An oscillation mode can be 
identified by the harmonic indices I and m (together with some statement about the radial eigenfunction) . 

For rotating stars, the problem is more complicated. First of all, rotation removes the degeneracy in m generating a 
richer spectrum than in a non-rotating star. Secondly, the identification of a mode with a single I is well defined only for 
non-rotating models. Nevertheless, each individual mode can be tracked by increasing the rotation rate step by step. It is also 
important to note that, even though the spectral properties of an oscillation mode change continuously along a sequence of 
rotating models, its overall parity is conserved. 

In this paper we discuss the three main classes of o scillation mode s. They are classified according to the dominant 
restoring force that acts on the perturbed fluid elements ( Cowlin j 1941). For polar modes the high frequency part of the 



spectrum, roughly above 1 kHz, is occupied by the fundamental f-mode and the acoustic p-modes. In stars where composition 
or thermal gradients are present, low-frequency g-modes are also present. These are restored by gravity through the buoyancy 
force. These various modes are non-trivial already in a non-rotating star. A polar mode labeled by the harmonic index I in a 
non-rotating conflguration is split by rotation into 21 + 1 distinct modes (identified by the value of m). 

In contrast, there are no non-trivial axial modes in a non-rotating fiuid star . The axial modes form a zero frequency 
subspace that describe stationary convective currents ( Lockitch fc Friedman 19991 ). In a rotating configuration, these modes 



become non-trivial thanks to the Coriolis force. Collectively, they are referred to as inertial modes. A general mode in this 
rotation-dominated class can have a mixture of axial and polar comp onents. However, since each mode has a distinct parity 
one can distinguish two separate classes of inertial modes. Following Lockitch fc Friedman ( 19991 ) we refer to these sets as 



polar- and axial-led. In the family of axial-led inertial modes, you find the so-called r-modes. They are particularly simple 
in that they are purely axial in the slow-rotation limit. For barotropic stars only the / = m r-modes exist, while in non- 
barotropic stars I 7^ m r-modes also appear. An interesting question concerns the fate of these r-modes in the barotropic limit. 
Some aspects of this question were resolved by Lockitch fc Friedman 1 19991 ). but interesting questions remain. Consider, for 



example, a weakly stratified star. In a fast spinning configuration the buoyancy should be dominated by the Coriolis force 
and the low frequency modes should essentially have an inertial nature. For lower rotation rates there should be a transition 
to a buoyancy dominated regime, where the g-modes become distinct. Of course, the relative strength of the Coriolis force 
and the buoyancy will vary throughout a star. As we will discuss in Section 17.31 one wou ld expect to find distinct r egions 
where buoyancy dominates the Coriolis force, and other regions where it is less important ( Dintrans fc RieutordI 200ol ) . It is 



then important to understand these two regions better. One would certainly want to know if it is possible to track individual 
modes from one regime to the other. Ideally, this would shed some light on the relationship between the inertial modes of a 
barotropic model and the g-modes of a non-barotropic system. 

Let us finally make a connection between the various classes of modes and our evolution problem. In the evolutions, we 
do not distinguish between axial and polar perturbations. Instead, we use the fact that perturbations of an axisymmetric 
background can be decomposed into two parity classes that satisfy different conditions in the equatorial plane. For the first 
class, which we will refer to as type I parity, SP^ , 5xp , fi^ , are all even under reflection with respect to the equatorial 
plane, while fg is odd. The I = m fundamental modes and the polar-led inertial modes belong to this class. Conversely, for 
the second class, that we will refer to as type II parity, 5P~ , Sxp , f,~, are odd and is even. Typical modes of a rotating 
star that belong to this parity class are the I + 1 = m fundamental modes and the axial-led inertial modes. 
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-1- R /R =0.65 




an In 

Figure 1. This figure shows the time evolution of the integrated flux component , defined in equation l|25| l. for two fast rotating 
models with Rp/ R^q = 0.75 and Rp/ Req = 0.65. Time is given in dimensionless units as a number of the star's rotational periods. These 
two simulations show the stability and the small numerical dissipation of our numerical code. 



6 THE CODE 

After decomposition with respect to 4>, the dynamical problem has two spatial dimensions, represented by the angle 6 and 
the new radial coordinate x = x {r,9). This radial coordinate is fitted to surfaces of constant pressure of the background star 
and it is normalised to unity on the equatorial radius. The choice of co ordinate simplifies the implementation of the surface 
boundary condition in a perturbed star considerably ( Jones et al. 2002h . The spatial variables {x^6) are then discretized on 



an evenly spaced two-dimensional grid, where ^ a; ^ 1 and ^ S ^ 7r/2. 

We evolve equations (|15|l - (|17| l in time with a Mac-Cormack algorithm, a second order accurate scheme that involves 
a predictor and a corrector step. In order to prevent spurious numerical oscillations we use a fourth-order Kreiss-Oliger 
dissipation with a parameter that depends on the numerical grid size. The perturbation equations are evolved for the internal 
points of the numerical grid, while boundary conditions are used to update the perturbations at the origin, the rotation 
axis, the surface and the equator. In order to use the fourth-order Kreiss-Oliger dissipation also near the grid boundary one 
ghost point is added to each grid slice. Numerical simulations are then carried out on [x, 6) grids with 60 x 32 points. Test 
simulations on finer grids with 60 x 64 or 120 x 32 points do not give significant improvements on the results. The numerical 
evolutions are stable for several hundred rotation periods. We extract the required mode frequencies from the evolutions by 
performing a Fast Fourier Transformation (FFT) of the time series of the perturbations. 

In order to have global information of the pulsations of the star, the perturbation variables can be integrated on the 
spatial grid (using a two dimensional integral). For instance, the integrated flux perturbation component is given by 

U= de f+rdr, (25) 

Jo Jo 

where R{0) is the radius of the star at any angle 9. This integrated quantity is shown in Fig[l] for two fast rotating stellar 
models with axis ratio 0.75 and 0.65, for evolutions lasting up to 120 rotation periods. This figure illustrates the codes stability 
as well as the very small numerical dissipation. Similar characteristics ar e observed fo r othe r perturbation quantities. These 
results demonstrate a significant improvement on the results reported in Ijones et al.1 |2002h . This improvement in stability 



allows us to analyse correspondingly longer stretches of data, providing more accurate localisation of modes in frequency 
space. 
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Q G f mode 1 = m = 2 

□ B f mode 1 = m = 4 



o o f mode 1 = m = 6 

^ ^ rmodel = m = 2 




n A I I < 1 < 1 1 1 < 1 < 1 1 1 1 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 

£l/(Gp^)"' 

Figure 2. This figure displays the frequencies of the I = m f- and r-modes as a function of the rotating rate for the Tp = 2 barotropic 
sequence of rotating stars. The mode frequencies are determined in the rotating frame. All frequencies are given in units of y/Gpc, where 
pc is the mass density at the stellar centre. T hese results test ou r numerical code in the limit of barotropic stars. The mode frequencies 
are in good agreement with those reported by I Jones et al. 



7 RESULTS 



7.1 Initial Data 



Once we fix the azimuthal harmonic index m the code evolves linear oscillations on an axisymmetric background. This means 
that modes with any value of I can be excited. This makes the problem different from a perturbation analysis in the frequency 
domain, where one tends to have some control over the specific angular dependence of the modes that are studied. The spectral 
content of the evolutions, which is extracted from an FFT of the time evolved perturbation variables, depends on the initial 
data used to excite the motion. 

On the initial time slice, a single oscillation mode can, in principle, be excited by providing the correct eigenfunctions for 
the perturbation variables. However, this is not a practical s trategy since it require s the mode-problem to be solved already. 
An alternative is to use an eigenfunction recycling method ( Stergioulas et al. 2004 : Dimmelmeier et al. 20061 ) ■ In this study 
we do not need these refinements. We are interested in a multi-mode analysis, so we can simply choose arbitrary initial 
perturbations that excite as many modes as possible in a single numerical simulation. 

We excite type I parity perturbations with a Gaussian radial profile in the pressure perturbation: 



5P = p exp 



ro 



{qR{0) 



Yn (0,^) 



(26) 



where R (9) is the stellar radius at polar angle 6, and the parameters ro and q determine the centre of the Gaussian and 
its width, respectively. The function Yn (9, 0) is the I — m spherical harmonic which describes a typical polar mode angular 
pattern of a spherical star. The initial conditions of the mass flux f and proton fraction Xp perturbations are set to zero. Type 
II parity perturbations are instead excited by the flux perturbation: 



p exp 



'■0 



(27) 



where Vjf (6, (fi) is a magnetic spherical harmonic ( Thornel IQSOl ). In this case, the pressure and the proton fraction perturba- 
tions both vanish on the initial time slice. Note that for both type I and type II parity, the coupling of the perturbed quantities 
in equations p5 p -(|17 p ensures these these simple choices of initial data sets do in fact excite non-trivial perturbations in all 
quantities. 



7.2 The Barotropic Case 

We first consider the problem for barotropic stars. This serves as a useful test-case for the reliability of our evolutions. The 
results also illustrate the improvements that we have made on the previous analysis ( Jones et al. 20021 ). 

In barotropic rotating stars the oscillation spectrum is formed by pressure, fundamental and inertial modes. Rotation 
removes the degeneracy (in m) of the non-axisymmetric pressure and f-modes and generates the new family of inertial modes. 
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Table 2. Inertial mode frequencies for barotropic rotating models with V p = 2 . The purely axial I = m = 2 r-mode is denoted by the 
symbol 'r and the inertial modes by 'i^, where k labels different modes w ith the same I. The second co lumn shows the overall mode parity, 
i.e. whether the mode is axial or polar-led. Frequencies determined bv lLockitch &: Friedmaiil ll 19991 ) are labeled with the acronym L-F 
and listed in the third column, while those extracted from our numerical code are shown in the last two columns for two slowly rotating 
models. In the formalism used by L-F, positive and negative cigenfrequencies correspond to counter and co-rotating modes respectively. 
In an FFT of a time evolved perturbation variable, they all appear positive, but their counter- or co-rotating character can easily be 
established by comparison with the L-F results and the dependence on the rotation of the star. 



Mode 


Parity 


L-F 


Rp/R^q = 0.983 


Rp/R^q = 0.95 


2r 


a 


0.667 


0.676 


0.695 


^ii 


P 


-0.557 


0.562 


0.549 


3i2 


P 


1.100 


1.101 


1.116 


*ii 


a 


-1.026 


1.021 


1.036 


*i2 


a 


0.517 


0.508 


0.536 




a 


1.358 


1.363 


1.375 




P 


-1.273 


1.291 


1.296 




P 


-0.275 


0.252 


0.258 




P 


0.863 


0.834 


0.884 




P 


1.519 


1.512 


1.519 


«il 


a 


0.422 


0.391 


0.409 




Figure 3. Rotating frame frequencies of polar-led inertial modes (left panel) and the r-mode and other axial-led inertial modes (right 
panel) as a function of the rotation rat e for V = 2 barotropic models. The inertial mode frequencies shown in this figure reproduce with 
good accuracy the results pr esented bv|jon es et al] |2002|) . In the limit of slow rotation, the mode frequencies are also in good agreement 
with the results reported bv lLockitch fc Friedm'anri|l999l) . see Table [2] 



which are restored by the Coriolis force. The main properties of the inertial modes are that they vanish linearly in the limit 
of zero stellar rotation and that they have mix ed nature. Even in the slow rotation limit the velocity perturbations of the 
inertial modes have both axial and polar parts ( Lockitch fc Friedman 19991 ). However, one can still identify two independent 
classes of inertial modes. The axial-led class with parity (—1)™+^ and the polar-led modes with parity (— 1)™. The velocity 
perturbation of an axial-led inertial mode can couple only with axial terms having / = m, m-\-2, . . . and with polar terms having 
Z = m -I- 1, m -I- 3, . . . . In the case of a polar-led inertial mode, the velocity perturbation couples only with I = m,m -\- 2, . . . 
polar terms and i = m-|-l,m-|-3, ... axial terms. A particular subset of the axial-led inertial modes is formed by the r-modes, 
which have purely axial velocity perturbations in the slow-rotation limit. For barotropic stars only I = m r-modes can exist. 
The non-axisymmetric modes of barotropic Newtonian rotating stars have been exten sively investigated in literature. We 



tested our numerical code by reproducing the results of the previous code developed by I Jones et al. 

I l|2002l) . finding a very 



good agreement. For stars with F/j = Ff = 2, we show in Fig.[2]the rotational splitting of the I = m f-modes and the I = m = 2 
r-mode. The final stellar model in the sequence is rota ting very fas t , at a bout 99 percent of the Kepler limit. The new code 
performs significantly better that the one discussed in I Jones et al.l (|2002l ). In particular, the improved stability allows us to 
study configurations that rotate closer to the break-up limit. Being able to carry out longer simulations, we can also improve 
on the precision of the extracted mode frequencies. 

In addition to the f-mode, many inertial modes are excited during a typical evolution. For type I parity perturbations 
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Table 3. Comparison of the first two I = 2,3,4 g-modes frequencies with the lYoshida Le3 l|200d ) results. The star is non-rotating 
and have adiabatic indices Tp = 2 and Ff = 2.0004. The frequencies are given in units of a j \/Gpc and the Yoshida and Lee values are 
denoted with YL. 



I 


gl 


gl 


g2 


g2 




YL 




YL 




2 


0.0189 


0.0187 


0.0129 


0.0123 


3 


0.0227 


0.0229 


0.0162 


0.0159 


4 


0.0255 


0.0252 


0.0189 


0.0188 



Table 4. The frequencies of the first three I = 2 g-modes for a non-rotating non-barotropic model with F^ = 2 and three different values 
of Ff = 2.05, 2.2, 2.4. The frequencies are given in units of cr/ yJGpc- This comparison of modes determined with a frequency (FD) and 
time domain (TD) approach provide a useful reliability check of our time evolutions. 





^gi 


^gi 


^g2 


^g2 


^g3 


^g3 




FD 


TD 


FD 


TD 


FD 


TD 


2.05 


0.195 


0.208 


0.134 


0.141 


0.102 


0.111 


2.2 


0.383 


0.396 


0.266 


0.271 


0.205 


0.209 


2.4 


0.531 


0.530 


0.370 


0.367 


0.287 


0.284 



and m = 2, we can identify polar-led I = 3 and I = 5 inertial modes. These are shown in the left panel of Fig. |3]for several 
barotropic rotating stars. For m = 2 type II parity initial conditions, axial-led I — 4 and I = 6 inertial modes and the I = m = 2 
r-mode are excited and their frequencies are shown in the right panel of Fig. [3l In orde r to ident i fy and test the inertial mode 
frequencies, we have compared them to the frequencies determined bv Ex)ckitch fc FriedmanI (199^. As is clear from the 
data in Table [2j the agreement is better than 6% for slowly rotating barotropic models with axis ratio Rp/Req = 0.983 and 
Rp/Req = 0.95. 

It is worth noting that, some of the eigenfrequencies in Table [2] are negative as these modes are co-rotating with the 
stellar fluid. Meanwhile, the frequencies that we extract will always be positive, as we are using an FFT of the time evolved 
perturbations. The co- and counter-r o tating character of the various modes can be established by comparison with the results 
determined bv lLockitch fc FriedmanI 1 19991 ). In some cases, the nature of the mode can also be inferred from the dependence 
on the rotation rate of the star. As can be seen from, for example. Figure [3] some of the co-rotating mode frequencies (in the 
rotating frame) tend to decrease at fast rotation rates. Meanwhile, no counter-rotating modes show this effect. 



7.3 Stratified Stars 

The oscillation spectrum of a stratified star is, in a sense, richer than that of a barotropic model. Fluid elements that have 
moved away from the equilibrium position can produce composition gradients in the stellar matter. For small displacements, 
buoyancy acts as restoring force by producing an oscillatory motion of the perturbed fiuid element. These oscillation modes are 
called composition gravity modes (g-modes). They must been discerned from thermal g-modes, which are due to temperature 
gradients. The latter are only relevant in newly born neutron stars, while the composition modes are relevant also in mature 
systems. Gravity modes have a polar nature, appear in the low-frequency band of the oscillation spectrum, typically around 
100 Hz, and exist even for non-rotating stellar models. Like f- and p-modes, each g-mode can be labeled by the harmonic 
indices I and m together with the number of nodes in the radial eigenfunction. We will label individual g-modes using the 
notation 'gk, where k labels different modes sharing the same value of the index I. Additionally, we will denote a retrograde 
(as opposed to prograde) mode as 'g]j 

There is an overlap in frequency space between the g-modes and the inertial modes. In fact, since both sets of modes 
emerge from the degenerate trivial modes of a non-rotating barotropic model, one might expect them to be intimately related. 
In order to understand this relationship better, we will investigate the low frequency oscillations of non-barotropic rotating 
neutron stars. The aim is to clarify whether inertial and g-modes are both present in the spectrum and to what extent it is 
possible to distinguish them. We also want to understand their behaviour in fast rotating models. 

These issues can be directly related to the dynamics of a perturbed fluid element and the balance between the two 
restoring forces that are acting on it, the Coriolis and the buoyancy force. In slowly rotating and stratified neutron stars with 
large composition gradients, buoyancy is expected to dominate the Coriolis force. As a result, a typical mode should appear 
as a g-mode. For higher rotation rates, the effect of the Coriolis force increases and may dominate the buoyancy force beyond 
a given stellar spin. In this regime, the oscillation mode should behave as an inertial mode. The threshold between these two 
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Figure 4. Results for non-barotropic rotating stellar models with = 2 and Ff = 2.05. The left panel shows the frequencies of the 
I = m = 2i- and g-modes, while the right panel shows the corresponding modes for Z = 3 and m = 2. In order to emphasize the rotational 
splitting of the g-modes and the onset of the CPS instability, the frequencies are shown in the reference frame of an inertial observer. 



regimes depends on the magnitude of the composition gradient and must obviously be smaller than the mass shedding limit 
of the star to be of any interest. 

The order of magnitude of the composit ion g-mode frequencies can be estimated by the Brunt- Vaisala frequency A^, 
which is defined by (see e.g., Unno et al. 19891 ) 

VP \ rv^ 



(28) 



For non-rotating stars this frequency can be written 



d\np 
dr 



(29) 



where the equilibrium configuration equations have been used to replace the pressure gradient with the gravitational acceler- 
ation. An averaged value of equation (|29p can be estimated by approximating the gravitational acceleration as (? ~ GM/R^ ~ 
(47r/3) GR{p), where (p) is the mass density average, and the density gradients as dp/dr ~ —pc/R. As a result, the averaged 
and dimensionless Brunt- Vaisala frequency N = N / yJGpc is given by: 

\ 1/2 

(30) 




We have tested equation (|30|) with the average of equation (|28p for the non-rotating and non-barotropic models used in this 
paper. The results agree to better than 5%. 

We want to study the effects of composition gradients on the spectrum in detail. We consider non-barotropic models 
with = 2 and frozen composition indices between Ff = 2.05 and Ff = 2.4 in order to investigate stars with low and high 
composition gradients. It is useful to begin by estimating at what rotation rate the Coriolis force will start to dominate the 
buoyancy. When Ff = 2.05, equation (|30|l suggests that A'' ~ 2f2 for a model with Q. / \JGpc ~ 0.16. For models rotating much 
faster than this, we would expect the Coriolis force to dominate. As a result, one would expect the low-frequency modes of 
such models to have a predominantly inertial character. As we will see later, the trend suggested by this rough estimate is 
brought out by the detailed mode-calculations. 

The non-barotropic evolution code has b een tested ag ainst two eigenvalue c odes developed by Yoshida fc Leel ( 200GI ) 
and Passamonti. Stavridis fc Kokkota3 ( 2008h . respe ctively. Yoshida fc Led (200oh solve for linear perturbations of slowly 
rotating neutron stars in Newtonian gravity, whereas Passanionti et alT ( 2008 ) study the pulsation spectra of slowly rotating 
relativistic stars in the Cowling approxi mation. We have comp ared the first two g-modes frequencies with the non-rotating 
and weakly stratified Newtonian model of lYoshida fc Led (2000|), which has F^ = 2 and Ff = 2.0004. The relevant frequencies, 
which are given in Table |3] show a very good agreement even though we are using the Cowling approximation. The second 
I — 2 g-mode have an error of four percent, while the others are accurate to within two percent. This difference is due the low 
frequency of the second 1 = 2 g-mode that requires a longer time evolution for extracting a more accurate v alue with an FFT 



We have further tested our code on neutron star models with larger composition gradients by using the lPassamonti et al 



(200: 



^ code. Since the time evolution code used in this paper is Newtonian, we have made a comparison between the two 



codes for a spherically symmetric background model with low compactness. For a polytropic model with F^ = 2 and a star 

-ll4 „-3 



mass Af = 1.17Mq and radius R = 16 km, giving a compactness ratio of M/R = 0.11. In terms of the dimensionless quantity 
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Figure 5. This figure displays the propagation diagram of two stratified models with T[ = 2.1 and Gpc = 0.084 in the left panel 
and Q/y'Gpc = 0.488 in the right panel. The Brunt- Vaisala frequency TV, the Lamb frequency Si and the Coriolis term 2flcos6 of the 
dispersion relation 13211 are given in yJGpc units. On the horizontal-axis the radial coordinate r is normalized with the stellar radius R. 
For simplicity, the N and Si frequencies are presented in the two panels for d = 7r/4, while three directions are shown for the Coriolis 
term, respectively 6 = 0,7r/4 and 7r/2. The two propagation diagrams illustrate the regions of the star where the various restoring forces 
dominate. Increasing the stellar spin, the Coriolis term becomes stronger than the Brunt- Vaisala frequency in almost all the volume of 
the star. Only in a small region near the equator, buoyancy can dominate even in rapidly rotating stars. 




o/(GpJ - 

Figure 6. This figure shows the Power Spectral Density (PSD) of the integrated flux perturbation component for a set of slowly 
rotating barotropic and non-barotropic models. All these models have axis ratio i?p/i?eq = 0.983 (corresponding to O/v'G'pc = 0.167), 
but different frozen composition adiabatic index Ff. The solid line corresponds to a barotropic star with Ff = 2. The counter-rotating 
^gi mode has been determined, from left to right in the figure, for stars with Ff = 2.05, 2.1, 2.2, 2.3, 2.4. This figure shows the dependence 
of the ■^gi mode frequency on the composition gradients. In the limit of Ff — > F^, the counter-rotating ■^gi mode seems to approach the 
•^12 inertial mode. This behaviour is clearer in Fig 



used in our time domain code, this model has Mj [pcR"^) — 1.262 which agrees very well (within one percent) with the non- 
rotating model in Table [T] By using these two approaches, we determined the g-mode frequencies of non-rotating stratified 
neutron stars in Tabled For non-barotropic stars with Ff — 2.05 the frequencies of the first three modes agree to better than 
eight percent, while for Ff = 2.4 they agree to better than one percent. This level of agreement is satisfactory given that the 
calculations are carried out within completely different frameworks. Thus we are confident that our time-evolutions are able 
to provide accurate g-mode results. 

The rotational splitting of the / = m = 2 f- and g-modes for stratified stars with Ff = 2.05 is shown, in the frame of an 
inertial observer, in Fig. |4] The left panel represents modes with type I parity, while the right panel shows results for type II 
parity. The two classes of modes clearly behave in a similar fashion. Results for other values of Ff are similar. Generally, we 
observe that the frequencies of the f- and g-modes tend to increase for stars with larger composition gradients, i.e. larger Ff. 
From the results in Fig. [4] it is easy to read off when the low-order g -modes become susceptible to the gravitational-radiation 
driven CFS instability I Friedman fc Schutj 19781 : Andersson 20031 ) . The instability condition is that the mode-frequency 
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/ ( G p J 



Figure 7. This figure displays the variation of the counter rotating ^gj mode with the adiabatic index Ff for a sequence of rotating 
stellar models. The polar-led inertial mode ^12 frequencies of the barotropic models with Ff = 2 are shown as a circle-solid line. This 
figure suggests that in rapidly rotating stars the counter rotating ^gj mode behaves similarly to the barotropic •^12 inertial mode. This is 
expected when the rotation rate is larger than the Brunt- Vaisala frequency A^, namely 2^2 > A'^, and the Coriolis force starts to dominate 
over the buoyancy force. 



0.12. This 



changes sign according to an inertial observer. The results in the figure suggest that this happens for Q,/^/G^ 
is in accordance with the typical values expected for the g-mode CFS instability. 

In fact, we can compare our result to an order of magnitude estimate based on the result of Lai ( 1999h . He estimates that 
the rotation rate of the g-mode instability onset corresponds to 

f7 = 0.68ao, (31) 
where o"o is the g-mode frequency for the non-rotating model. From this we find that one would expect the instability of the 
^gi mode for the Ff — 2.05 model to be triggered in a star rotating at ^Gpc ~ 0.14. This is close to the (presumably more 
accurate) value determined with our numerical code. 

However, the g-mode instability is not thought to be very important since the g-modes are not radiating gravitational 
waves efficiently. Viscosity is expected to suppress the instability in more realistic models. In comparison, the f-mode, which 
radiates gravitational waves efficiently does not go unstable until the star spins extremely fast. For stellar models with 
F^ = 2, the Z = m = 2 f-mo de is not expected to become unstable before the mass shedding limit rotational rate is 
reached ( Ipser fc Lindblom 1990l ). As is clear from Fig. [l] the f-mode instability point occurs for faster spins than we are able 
to study. This accords well with previous results in the literature. 

Let us now study the low frequency band of rotating non-barotropic stars. A local analysi s of uniformly rot ating stars in 
the Cowling approximation leads to the following dispersion relation for low-frequency waves (|Unno et al.lll989l ): 



iv^fci -H (2n ■ k)^ 

fe2 



(32) 



where k is the wave vector and fc^ is its comp onent orthogonal to the apparent gravity. These modes are generally referred to 
as inertia-gravity waves bv lUnno et Zl (Il989l ). as they are restored by both gravity and the Coriolis force. In the barotropic 
case, N = 0, equation (|32p describes the inertial waves, while for a non-rotating model it reduces to the usual gravity waves. 
In stratified and rotating neutron stars, it is thus to be expected that the low-frequency oscillation modes have an hybrid 
character with gravity dominating at low rotation rates and the Coriolis force taking over at fast stellar spin. 

However, it is worth noting that the buoyancy and the Coriolis force can be locally dominant in different regions of the star. 
This can be understood from a propagation diagram that illustrates the spectral properties of the star. For instance, we consider 
two stellar models with frozen adiabatic index Ff = 2.1 and with angular velocity Q,/^Gpc = 0.084 and Q./\/Gpc = 



respectively. For these models we determine the Brunt- Vaisala frequency A'^ (|28p . the Lamb frequency ljUnno et al 



198 



i) 



s = 



[0+1)^2 



(33) 

and effect of the Coriolis force on the dispersion relation (|32p . i.e. 2ficos6'. The role of the Lamb frequency is to estimate the 
propagation region for acoustic modes. For these two stellar models, the propagation diagram for the I = m = 2 modes is 
shown in Fig. [5]using the normalized radial coordinate r / R. For the slowly rotating star with Rp/Req = 0.996, the N and Si 
frequencies do not vary much, while the effect of the Coriolis force term 2S7cos6' changes according to its angular dependence, 
from zero at the equator to the maximum at the pole. In this case, the star is spinning so slowly that the buoyance (N) 
is dominant in most of the star. Only near the centre, for r <; 0.2R, is the Coriolis force stronger. As a result, we should 
expect the inertia-gravity modes to behave as g-modes. Later, this will be confirmed by a global mode analysis. When the 
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Figure 8. Eigenfunctions of the flux component perturbation for rotating stellar models with rotation rate Q/y/Gpc = 0.084 (left 
panel) and fl/^/Gpc = 0.488 (right panel). The yellow surfaces refer to the '^12 inertial mode of barotropic models with Ff = Fp = 2, 
while the red surfaces represent the counter-rotating ^gj mode of stratified stars with Ff = 2.1 and F^ = 2. This figure illustrates how 
g- and inertial mode eigenfunctions become similar for increasing rotation rates. 
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Figure 9. Eigenfunctions of the flux component perturbation for the rotating stellar model with rotation rate fl/y/Gpc = 0.084. 
The left panel shows the contour plot of the ^12 inertial mode for the barotropic model with Ff = F^ = 2. The right panel displays the 
contour plot of the counter-rotating ■^gi mode of a stratified star with Ff = 2.1 and F^ = 2. These plots give a complementary view of 
the eigenfunctions shown in the left panel of Fig. ((Sj, and confirm the difference between the two modes in a slowly rotating model. 
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Figure 10. Eigenfunctions of the fiux component perturbation for the rotating stellar model with rotation rate Q/^/Gpc = 0.488 
The left panel shows the contour plot of the '^12 inertial mode for the barotropic model with Ff = F^ = 2. The right panel displays the 
contour plot of the counter-rotating ^gi mode of a stratified star with Ff = 2.1 and F^ = 2. These plots confirm the results of Fig. ^ 
showing that the structure of these two eigenfunctions becomes similar in fast rotating models. 



14 A. Passamonti et al 




Figure 11. Frequencies of some selected g-modes (solid lines) for Ff = 2.05 models and polar-led inertial modes (dashed lines) of 
the barotropic F^ = Ff = 2 model. All frequencies are measured in the rotating frame. The dominance of the Coriolis force in rapidly 
rotating stars is evident from this figure. We see that the g-modes are mainly restored by the Coriolis force when H/^/Gpc > 0.3, and 
thus become similar to the barotropic inertial modes. 



star rotates faster, as in the model with Rp/Req — 0.85 (right panel of Fig.[5ll, the Coriolis term is more important and larger 
than the Brunt- Vaisala frequency in the most of the star. In order to keep the figure clear, we have shown the and Si 
frequencies only for 9 — 7r/4. The change in the other directions is visible, but does not modify the qualitative properties 
of the propagation diagram. However, even in weakly stratified and rapidly rotating stars there is a region near the equator 
where buoyancy dominates over the Coriolis force. This region, which is due to the cos^ dependence of the Coriolis term, 
shrinks for higher rotation rates. 

The inertia-gravity modes can be identified in our time domain simulations by studying th e FFT of the pe r turba - 



tion evolution and extr a cting the ID and 2D mode eigenfunctions with the code developed by IStergioulas et al.l (|2004l ) 



and iDimmelmeier et al.l (|200q ). In order to study the dependence of the modes on the stellar spin and stratification, we 



construct a set of equilibrium models by gradually varying the two parameters and Ft. We can then follow the change of the 
mode frequency and eigenfunction shape and correctly distinguish it among the other peaks in the spectrum. Let us consider, 
for instance, a sequence of slowly rotating models with fixed rotation rate Q,/y^Gpc = 0.167 and different Ff. Typical FFT 
results are shown in Fig. [B] where we have identified the ^i2 inertial mode for the barotropic model and the counter-rotating 
■^gi mode for a set of stars with increasing composition gradients. The results suggest that the g-mode frequency approaches 
the barotropic inertial mode in the limit Ff — > 2. This association becomes more evident if we combine these results for a 
set of rotating models, see Fig. [7] When the star rotates with 20, > A'^, the retrograde ^gi mode seems to behave as the 
barotropic ^12 inertial mode. This transition obviously depends on the level of stratification. This property is also noted in the 



mode eigenfunctions. In Figs. [8 H10I we show how the eigenfunctions of the mass fiux component of the counter-rotating 
mode approach the shape of the "^12 mode when 20, > A'^. In this case, the non-barotropic star has Ff = 2.1 and an 
estimated Brunt- Vaisala frequency (N/^Gpc) — 0.447, see equation (|30p . In Fig. [8] we provide a three-dimensional picture of 
the eigenfunctions, for rotation rate Q/^Gpc = 0.084 in the left panel, and for Oj^JGpc = 0.488 in the right panel (the latter 
corresponding to a star with an angular velocity slightly larger than the Brunt- Vaisala frequency). Clearly, for the slower 
rotation rate (left panel) the the ^gi eigenfunction is quite different from the eigenfunction of the barotropic "^12 mode. In 
particular, the nodal curve of the two modes do not coincide near the equator. In contrast for the more rapidly rotating star 
(right panel) the nodal curves of the two modes are very similar. This can be seen particularly clearly in the contour plots of 
the eigenfunctions given in Fig. H (for O/y/GfTa = 0.084) and Fig. [10] (for O/^/Gp'a = 0.488). 

When we extend this approach to the other g and inertial modes, we find further mode associations. Combined results 
are provided in Figs. 1111 and [121 These figures display the spectrum of barotropic models and stratified stars with Ff = 2.05. 
The results show that individual non-barotropic g-modes tend towards specific barotropic inertial modes as the rotation rate 
of the star is increased. This agrees with the notion that the Coriolis force dominates the buoyancy in fast spinning systems. 
From the results in the two figures we can identify various g-modes with barotropic inertial-mode counterparts. For type I 
parity and m = 2 this leads to: 

4 r 5. 4 p 5. 2 r 3. 2 p 3. 2 r 5. 2 p 5. /o/l^ 

gl < > 14, gl< > ll, gl < > 12, gi < > ll, g2 < > 13, g2 < > 12, (34) 

where the upper index r and p denotes retrograde and prograde g-modes respectively. Meanwhile, for type II parity and m — 2 
we find that: 
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Figure 12. Frequencies of some selected m = 2 g-modes (solid lines) for Ff = 2.05 models and axial-led inertial modes (dashed line) of 
the barotropic models F^g = Ff = 2. Frequencies are measured in the rotating frame. The counter- and co-rotating ^gi modes approach 
the frequencies of the axial-led inertial modes in the limit of fast rotation. In contrast, the r-mode is essentially unaffected by composition 
gradients. 



gl < — > 13 , gi < — > u • (35) 

From Fig. [12] we also see that it is only the I — m r-mode that does not liave a g-mode counterpart. As expected, the r-mode 
is essentially the same in the barotropic and non-barotropic models. In the type II parity perturbations we could accurately 
determine only the first g-mode, as the r-mode pe ak was dominant in th e FFT and overwhelmed the higher order g-modes. 

These results are consistent with the study of lYoshida &: Lee where the authors worked within the slow rotation 

approximation and in the frequency domain. They studied inertia-gravity modes of F/j = 2 polytropic models with weak 
stratification and slow rotation. In fact, their frozen adiabatic index corresponds to Ff = 2.0004 and the fastest rotating 
model they consider represents Q,/ [GM / — 0.1, where M and R are respectively the mass and radius of the non- 
rotating models. In our units this rotation rate in equivalent to Q,/{GpcY^^ = 0.113 and an axis ratio Rp/Req > 0.992. 



In contrast, in the present work we investigate rapidly rotating models up to i}/(Gp. 



1/2 



0.717 {Rp/Reg = 0.6) and with 



larger composition gradients. We have tested our results against those of Yoshida fc Leel ( 2000l ) by numerically evolving m = 2 
perturbations of the two slowest rotating models reported in Table [1] with Ff = 2.0004. The fre quencies have an accu racy to 
better than 2 percent. Finally, equations (|34|) and (|35|l respect the connection rules reported in Yoshida fc Leel 1 2000l ). 



8 CONCLUDING REMARKS 

Information concerning the oscillation spectra of realistic neutron star models helps develop our understanding of the physics 
required to describe these objects. Any spectral property can in principle be attributed to particular physical quantities or 
configurations of the star. Astero-seismology studies, using either electromagnetic or gravitational signals (or indeed both) 
may thus help constrain the state of matter at supernuclear densities. Of course, in order to facilitate this, we need to improve 
our theoretical models and clarify the origin of various stellar pulsation features. 

In this paper, we have studied the pulsations of stratified and rapidly rotating neutron stars with the aim of understanding 
the dependence of the composition g-modes on the rotation rate of the star. The stellar pulsations were studied using the 
linearised Euler and conservation equations on an axisymmetric background. In order to simplify the problem, we used the 
Cowling approximation where the gravitational potential perturbations are neglected. This approximation generates only a 
small error in the g-mode frequencies. The system of perturbation equations was evolved in time by a 2D numerical code based 
on a standard finite differencing scheme. This code was tested against various results available in literature and demonstrated 
good accuracy. 

Since both Coriolis and buoyancy forces act on a perturbed fiuid element of a rotating and stratified star, the low- 
frequency modes have a mixed inertia-gravity character. They typically behave as g-modes in the slowly rotating limit, while 
they assume the properties of the inertial modes when the star rotates rapidly. By comparing the oscillation frequencies and 
the associated eigenfunctions of barotropic inertial modes and the inertia-gravity modes of non-barotropic stars, we have 
demonstrated how the two sets of modes are associated at fast rotation rates. 

Our results show that it would be difficult to deduce the presence of composition gradients from the inertia-gravity mode 
spectrum of fast rotating stars. However, as the neutron star ages and spins down, the dependence of the inertia-gravity modes 
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on the star's rotation and their deviation from the inertia! barotropic modes can, at least in principle, be used to estimate 
the Brunt- Vaisala frequency and the degree of stratification. 

As an extension of this work, we are currently developing a numerical code for studying the dynamics of rapidly rotating 
superfluid neutron stars. This is an important development since all mature neutron stars are expected to have superfluid 
components i n the core. I t is al so well known that the associated multi-fluid dynamics leaves an imprint on the stellar oscillation 
spectrum, see lLin et al.l l|2007l ) for a recent discussion. We hope to be able to report on the initial results of this investigation 
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